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Abstract. We present a simple method to solve spherical harmonics moment systems, 
such as the the time-dependent Pjv and SPn equations, of radiative transfer. The method, 
which works for arbitrary moment order TV, makes use of the specific coupling between the 
moments in the Pjv equations. This coupling naturally induces staggered grids in space and 
time, which in turn give rise to a canonical, second-order accurate finite difference scheme. 
While the scheme does not possess TVD or realizability limiters, its simplicity allows for 
a very efficient implementation in Matlab. We present several test cases, some of which 
demonstrate that the code solves problems with ten million degrees of freedom in space, 
angle, and time within a few seconds. The code for the numerical scheme, called StaRMAP 
(Staggered grid Radiation Moment Approximation), along with files for all presented test 
cases, can be downloaded so that all results can be reproduced by the reader. 



1. Introduction 

The purpose of this paper is to present a simple, yet accurate solution method for the P/v 
equations of radiative transfer, and its efficient implementation in Matlab. The key idea is 
to make use of the specific coupling of unknowns that is induced by the spherical harmonics 
being a family of orthogonal polynomials. This leads to a natural staggered grid on which 
the equations are discretized. 

The Pjsr method (cf. [lj) is one of several ways to discretize the equation of radiative 
transfer. It is often introduced as an approximate method (method of moments) to reduce 
the high dimensionality when the kinetic equation of radiative transfer, which is formulated 
on a six-dimensional domain (one time, two angle, three space), is discretized. Another way 
of interpreting the P/v equations is to view them as a spectral discretization in the angular 
variable. 

The efficient numerical solution of the Pat equations has become a recent subject of interest 
[HI QjjJ Q31 Q] . The Pjv equations have several advantages over other more direct discretiza- 
tions, such as discrete ordinates, most prominently rotational invariance. The lack of this 
property leads to the so-called ray effect in discrete ordinates approximations (cf. [H]). The 
key property that the numerical method presented in this work is based upon, is also exclusive 
to spherical harmonics moment methods, namely a specific coupling structure between the 
moments. The main drawback of the P/v equations is that they, being a spectral method, 
can exhibit Gibbs phenomena, i.e., oscillatory behavior that is not present in the solution 
of the original kinetic equation. In the context of radiative transfer, this can yield negative 
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and therefore unphysical particle densities. In many cases, the Gibbs phenomenon is not a 
major problem (given the oscillations are small in amplitude), since the P/v equations remain 
always well-defined. However, sometimes the unphysical particle densities are unacceptable, 
and several recent works have addressed this fact [7J Q3] . 

The numerical method presented here does not possess limiters in the hyperbolic solver 
nor does it overcome the Gibbs phenomenon. However, its simplicity allows for fast and very 
highly resolved computations, so that one can usually reduce the spurious oscillations to an 
acceptable magnitude by choosing the moment order N sufficiently large. 

This paper is organized as follows. In Sect. [2] we introduce the P/v equations in slab 
geometry and point out their structure. This is done for didactical reasons, because the 
notation for the P/v equations in two space dimensions (derived in Sect. [3]) becomes quite 
tedious. The staggered grid method is presented and analyzed in Sect. |4j and in Sect. [5] the 
efficient implementation of the numerical scheme in Matlab is outlined. Numerical examples 
are presented in Sect. [6j We attempt to meet the standards of reproducible research in the 
computational sciences, laid out by LeVeque [11]. The source code of our package StaRMAP 
(Staggered grid Radiation Moment Approximation), along with files to generate all this paper's 
figures, as well as additional examples, are available to the reader online [21]. 



2. The Slab Geometry P/v Equations 
We consider the radiative transfer equation in the form [2] 

d t ip{t, x,Sl) + Sl- V x ip(t, x, SI) + E t (t, x)ip(t, x, SI) 

f (!) 
= / £ s (i, x, SI ■ SI )ip(t, x, SI ) dSl + q(t, x, SI) . 

Js 2 

The quantity tp, which is defined for time t > 0, space coordinate x G M 3 , and direction 
SI £ S 2 , is the density of photons that undergo scattering and absorption in a medium. The 
medium is characterized by the absorption cross section S a , the scattering kernel S s and the 
total cross section = T, s q + S a (where £ s o is defined below). In addition, there is a source 

q- 

The slab geometry radiative transfer equation is obtained by considering a slab between 
two infinite parallel plates. Assume for instance that the z-axis is perpendicular to the plates. 
If the setting is invariant under translations perpendicular to and rotations around the z-axis, 
then the unknown tp depends only on the z-component of the spatial variable, and one angular 
variable [i (cosine of the angle between direction and z-axis). 

To obtain the P/v equations, we express the angular dependence of the distribution function 
in terms of a Fourier series, 

oo 

iP(t,z, fJ ,) = J2Mt,z)^Pi(») , (2) 

£=0 

where Pi are the Legendre polynomials. These form an orthogonal basis of the space of 
polynomials with respect to the standard scalar product on [—1, 1]. 
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One can obtain equations for the Fourier coefficients 



ipt 



(3) 



by testing the radiative transfer equation ([!]) with Pg and then integrating. Thus we obtain 
(suppressing the arguments) 



d t ip£ + d z / fxP£ip d/u + Ytti^i 



for I = 0, 1, . . . , where 

^« = S 4 — £ s 



sO 



and £ 



2?r 



P £ (fj,)E s (fj,) d/i . 



-i 



Two properties of the spherical harmonics are crucial for our method. First, we observe that 
by the procedure above we have diagonalized the scattering operator on the right hand side 
(the Legendre polynomials are eigenfunctions of the scattering operator). Second, a general 
property of orthogonal polynomials is that they satisfy a recursion relation. For instance, the 
Legendre polynomials Pe satisfy 

^Pe(n) = 2£3^-i(M) + im p i+^) ■ 



Using this fact and truncating the expansion at 
equations 



N we arrive at the slab-geometry Pj 



N 



d t ip£ + d. 
This system can be written as 



(4) 



dtu + M ■ d z u + C ■ u = q , 



where 



M 



/0 1 



v 



2 
3 


3 
7 



and C 



7 



f0 



\ 



The two properties mentioned above lead to 

Lemma 1. The time derivative of ipe for even (odd) t depends only on the spatial derivative 
ofipk for odd (even) k, and on the value of ipe itself. 

Lemma [T] creates an analogy to the wave equation (with decay), and thus motivates a 
discretization of the slab geometry P/v equation Q on staggered grids, i.e., all the components 
with odd i are placed in the middle between the components with even £, and the spatial 
derivative is approximated by central differences. The numerical scheme presented in Sect. [4] 
generalizes this analogy in the two-dimensional case. 
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3. The Two-Dimensional P n Equations 



In this section, we adopt the notation and the form of the Pn equations as in [Tj. The 
complex- valued spherical harmonics are defined as 



Yr(f,,d>) = (-iy 



21+1 (<-m)l jmtp™^ ; 



4tt (e+m)\ 

where £ > and — I < m < I. Here, Pf 1 are the associated Legendre polynomials. The 
spherical harmonics form an orthonormal family on the unit sphere. They satisfy a recursion 
relation of the form 



m—l v m-l | jm-lym-1 , m+l-i^m+1 



1 l-l 



€+1 £fl 



jm-lym-1 , m+l-i^m+1 

u e+i 1 e+i 1 1-\ 



0( n m ~V m _1_ A m V m "1 

*\ a e-i I e-i u e+i I e+i) 



with the coefficients [1] 



(l-m+l)(l+m+l) 
(2^+3) (2^+1) 

(£-m)(£-m-l) 
(2£+l)(2f-l) ' 



(£-m)(^+m) 
(2£+l)(2^-l) ' 

(l-m+l)(l-ro+2) 
(2<?+3)(2£+l) 



f r, 



•^+1 I e+i 



(l+m+l)(l+m+2) 
(2^+3)(2£+l) 

(e+m)(e+m-l) 
(2£+l)(2£-l) • 



This form already shows a pattern in the coupling of the different moments, that is similar 
to the slab geometry case. 



We multiply by Y™, integrate over f2, and define the expansion coefficients 



As in the slab geometry case, the scattering term becomes diagonal 



5-' 



Y £ m (n) / • n')ip(t, x, n') an' do = s a ^(t, x) . 



s 2 



Altogether, we obtain the well-known complex-valued Pn equations 

+ 5 z (a, m _ 1 ^i + &m)+S«^ 
for < ^ < oo and —£< m < £. 



,171 — 1 „ I.TYl — 1 
1 



jm-1 / m-1 , 



m+1 _ ftn+1 im+l\ 

i v^-i -/m-i Ve+i ) 



(5) 



In this work we consider the two-dimensional real-valued Pn equations, which we now 
derive. There is, however, no conceptual difference to the three-dimensional equations. The 
reduction is again done via symmetry. This means that we actually solve three-dimensional 
radiative transfer, but in a geometry that reduces the number of unknowns. For a two- 
dimensional domain D, consider the infinite cylinder D x R C M 3 . We take the angular 
variable to be aligned with the z-direction, 

= ( V 1 — yu 2 cos 4>, \/l — fi 2 sin (j), /u) T . 

If we assume that all data (coefficients, initial and boundary conditions) are ^-independent, 
then the solution tp is z-independent and additionally an even function in /i. Therefore, if 
£ + m is odd, the associated Legendre polynomial Pl a is an odd function in /x, and as a 
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consequence the moments for which £ + m is odd have to vanish. Thus we are left with the 
(still complex) moments 

We thus obtain the following matrix formulation of the P/v equations 

^r 1 -fl 

















+ d x \ 


^2 









d 



-fS 



-fl 



o 



+ dy\ 



-fl 



-ci 



-f°2 



-fl 



r 1 
c l 



* 
* 
* 

>o ' 

^2 
Y>1 



>0°' 

v> 2 " 2 
^ 2 



+ 



->ti 



H2 



->t2 



£* 2 





r -08 1 




r«8i 
























^ 2 








q°2 











We call the matrix behind the x-derivative (including the \) M£ omplex , and respectively the 
matrix behind the y-derivative (including the |) My° mplex . We denote the matrix containing 
the T, t £ by C. 

The last step is to transform this system to real variables. Note that 

= (-i)"*Vf • 

Real variables R™ (for < m < I) and I™ (for < m < t) can be obtained by setting 

fl? = V? 

and for m / 

(-1)" 



rm 



/™ = t^(^r - (-i) m v7 m ) • 
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The factor (— l) m is chosen so that the coefficients can be compared more easily to well-known 
moments. For example, 

rl r2n 

m)dn = ^= / y(/^)#d/i, 



Rq 



1 

/An 



An 



1 JO 



R\ 



P , /*1 f''2lY 

/ n x ip{n) dn = J-^ \ \ V^-v 2 cos 0^,0) , 

J An V J-l JO 



J4vr J-l JO 

The factor \/2 is chosen to make the transformation from ip™ to R™, l™ unitary. If we encode 
this linear relationship into a matrix S, so that 













R\ 








I{ 






u : = 


R 2 
T 2 
R 2 


= S- 


TP°2 



then we obtain the real- valued system matrices for the P/v equations as 
MI cal = SM^P^S" 1 and M rcal = SMr^^" 1 



(6) 



For example, the P3 matrices are 

V2d^ 





M rcal = 1 
lvl x 2 



and 



M. 



real 1 





V2d^ 



d 



-2 
2 







d 2 2 



-\/2/ 2 ° 



Vld^ 1 



l 2 





-d 2 2 





d 2 2 










d 2 2 









df 
















d 3 3 
















-/3- 1 





V2d^ 





-/3- 1 








d 2 2 





d- 2 2 











-d 3 - 3 





df 











-/3- 1 





/3- 1 





y/2d^ 



-fz 1 

fg 3 

V2d^ 






df 





d 3 3 



















-/3- 1 





/3- 1 



Vld^ 1 
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From the previous calculation it can be seen that the matrices M£ eal and ilPj eal couple the 
variables RJ 1 and I™ in a specific way, as follows. 

Lemma 2. The following variables are coupled, provided that they are defined (i.e., for R™ 
we must have £ >0, < m < £, for IY 1 we must have £ > 0, < m < £): 

• The time-derivative of R™ depends only on the x-derivatives of Rl^-x , on the y- 
derivatives o/PZ^ 1 , and on R™ itself. 

• The time- derivative of I™ depends only on the x-derivatives of PZ^ 1 , on the y- 
derivatives of R^^ 1 , and on I™ itself. 

The scattering matrix C is diagonal and block-wise constant, and therefore invariant under 
the transformation Consequently, the real- valued P/y equations read 

d t u + M v x eal ■ d x u + M T y cal -dyU + C -u = S -q, (7) 

where S ■ q contains the real-valued moments of the source q. We also note that M£ eal and 
M^ eal are both symmetric. 

Remark 1. The simplified P/v (SPn) equations derived in [T7] have the same coupling pattern 
and can thus be solved with the same numerical scheme, presented in Sect. |4j as the P/v 
equations. 



4. Numerical Method 

We now develop a numerical scheme for linear systems of hyperbolic balance laws of the 
form 

d t u + M x ■ d x u + M y -d y u + C-u = q, (8) 

where the matrices M x , M y , and C possess very specific patterns of their nonzero entries 
that admit the systematic placement of the components of the solution vector u(x,y,t) on 
staggered grids. Specifically, let the solution of Q have M components, i.e., u(x,y,t) £ IR , 
and the source q(x,y,t) G R M and the matrices M x ,M y G R ArxM and C(x,y,t) G R MxM are 
of appropriate sizes. Moreover, the matrix C(x,y,t) is diagonal, and the matrices M x and 
My are constant-coefficient, and possess patterns of their nonzero entries, as described in 
Lemma [2| Hence, the real- valued P/v equations ([7]) are covered, as are the SPn equations 

m 

4.1. Spatial Approximation on Staggered Grids. We consider the partial differential 
equation Q to hold in the interior of a rectangular computational domain Q = (0, L x ) x 
(0, L y ), and on each of the two boundary directions (horizontal and vertical), we allow for 
one of two types of boundary conditions: periodic or extrapolation, as described in more 
detail below. The domain 0, is divided into n x x n y rectangular cells of size Ax x Ay, where 
Ax = L x /n x and Ay = L y /n y . The center points of these cells then lie on the grid 

Gil ={((*- l)&x, (j - \)Ay) | i € {1, ... , n x }, j G {1, . . . , n y }} . (9) 

We always place the first component of the solution vector, i.e. the scalar flux Pj], on this 
cell-centered grid G\\. As an example, Figure [l] shows the division of the rectangular domain 
(light gray) into a 5 x 3 arrangement of cells. The grid G\\ is depicted by gray circles. The 
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Figure 1. Staggered grid of 5 x 3 grid cells, with periodic b.c. in the x-direction, and extrapolation 
b.c. in the y-direction. Shown are solution grid points (black boundaries), periodic extension points 
(light blue), and extrapolation ghost points (light red). 



key principle of the numerical scheme is to place the remaining solution components on grids 
that are staggered with Gn. 

To reiterate, the condition on the nonzero entry patterns of M x , M y , and C, given in 
Lemma [2j can be reformulated as follows: the components of u can be distributed into 
four disjoint sets, according to {1,2,..., A/} = In U/21 U J12 U J22, such that the following 
properties hold: 

(M x ) itj = V(t,i) i {{hi x hi) U {hi x In) U (J12 x I22) U (I22 x J 12 )) , 

{M y )ij = V(i,j) i {{hi x I12) U (J12 x In) U (J 2 i x / 22 ) U (J22 x hi)) , (10) 

dj = V(i,j) £ ((7ii x In) U (I 2 i x / 2 i) U (J12 x J 12 ) U (J 22 x I 22 )) • 

With this distribution of the indices of the solution components, we consider four fully stag- 
gered grids: Gn, defined above, and in addition 

G 2 i = {{iAx, {j - |) Ay) I i G . . . , nj, j G {1, . . . , n^}} , 

Gi 2 = {((i - |) Ax, j'Ay) | i G {1, . . . ,n x }, j G {p y ,.. . ,n y }} , (11) 

G 22 = {(iAx, jAy) I j G {fz, . . . ,n z }, j G {pj,, . . . ,n y }} , 



where 



Pa 



(12) 



!0 extrapolation b.c. m x I extrapolation b.c. m y 

and p y = < 
1 periodic b.c. in x 11 periodic b.c. in y 

In Fig. [TJ the grid G21 is depicted by gray top-pointing triangles, the grid G12 by gray right- 
pointing triangles, and the grid C?2 2 by gray squares. 

Having defined these fully staggered grids, the solution components with indices in he 
are assigned to the corresponding grid Gki, where k,£ G {1,2}. On these staggered grids, 
spatial derivatives of a function w can be approximated by the half-grid central difference 
approximations 

d x w{iAx,jAy) « ^ {w{{i + I) Ax, jAy) - w{{i - \)Ax,jAy)) V*,i € |Z , 
d y w{iAx,jAy) « ^ (w(»Aa:, (j + \)Ay) - w{iAx, {j - \)Ay)) G \TL , 



(13) 
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and we denote the resulting finite difference operators D x and D y , respectively. With these 
finite difference approximations, ^-derivatives of components on the grid G k i are associated 
with the grid G^-k,i, and y-derivatives of components on the grid G k i are placed on the grid 
Gk3—i, where k,£ G {1,2}. The nonzero entry patterns (10) guarantee that the distribution 
of the indices of u into the sets In, I21, I12, and I22 is exactly reproduced by the distribution 
of the indices of M x ■ D x u + M y ■ D y u + C ■ u. Moreover, the components of the source vector 
q in (Tsl) are placed on the same grids as the corresponding components of the solution vector. 



4.2. Boundary Conditions. In each of the two coordinate directions, we prescribe one 
of two types of boundary conditions. In the x-axis direction these conditions, and their 
implementation, are as follows: 

• Periodic: The solution satisfies u(x + L x ,y,t) = u(x,y,t). On the staggered grids, 
periodicity is implemented by "wrapping around" the grid data, by defining: for 
all components k G G21 U G22 set Ufc(0, jAy) := u k (L x , jAy) Vj G I y , and for all 
components k G G\\ U G\2 set Uk(L x + | Ax, jAy) := u k (^Ax, jAy) Vj G I y , where 

Iy = i 2,Py ' 2~Py + 2 > ■ • • > n y ~ 2 > n y } • 

• Extrapolation: All solution components with k G G\\ U G\2 satisfy d x Uk = at the 
left and the right boundary. On the staggered grids, this condition is implemented by 
constant extrapolation, i.e., by "copying" grid data onto ghost points, by defining: 
u k (-\Ax,jAy) := u k {\Ax,jAy) Vj G I y , and u k {L x + \Ax,jAy) := u k (L x - 
\Ax,jAy) Vj G I y , where I y = {\p y , \p y + \, . . . ,n y - \,n y }. 

The definition and implementation of boundary conditions in the y-axis direction works anal- 
ogously. Figure [l] shows an example with periodic b.c. in the x-direction (grid data that is 
wrapped around is depicted in blue) and extrapolation b.c. in the y-direction (ghost points 
are depicted in red). 



4.3. Time Stepping. The time-derivative in ^ is resolved by bootstrapping, i.e., data 
on the grids G\\ U G22 is updated alternatingly with data on the grids G21 U G\2- This 
approach is natural, since in the approximate advective part of Q, i.e., M x ■ D x u + M y ■ D y u, 
the solution components in In U J22 (called the "even components") depend solely on the 
components in I21 U I12 (called the "odd components"), and vice versa. One possibility to 
perform bootstrapping is by staggering the data in time, i.e., when using a time step At, the 
even components would be defined at times nAt, while the odd components would be defined 
at times (n + |)Ai. While such a space-time-staggered approach yields very elegant update 
rules, the placement of different solution components at different times causes inconveniences 
in the implementation of initial conditions, as well as the evaluation of the full solution at a 
fixed time. Here, we circumvent these problems by defining all solution components on the 
same time- grid nAt, and conducing a time step from t to t + At via a Strang splitting |22j of 
sub-steps, as follows. 

For the purpose of the presentation in this section, consider the solution vector be written 
in block form 



where n 6 is the vector of the even components, and u° is the vector of the odd components. 
Moreover, let the same block form apply to the source vector q and the matrices M x , M y , 
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and C. Then, with the central difference approximations (13), equation ([8]) turns into 

[tpi r o m-1 r^i r o m-i r^i rc e o 1 T^i rri n41 

' [u°\ ^ [M° c J x [v°\ ^ [M° c J y [u°\ ^ [ C°\ ' [u°\ [q°\ ' { 1 

We now define two evolution operators: one that updates the even components only, while 
"freezing" the odd components; and one that updates the odd components only, while "freez- 
ing" the even components. To advance the even (odd) components from t to t + At, we 
consider the odd (even) components to be constant on the interval [t,t+ At]. The same holds 
true for the matrix C and the source vector q, both of which are evaluated at time t + \At. 
With these assumptions, the block system (14) decouples into two ODEs 

{ dtif + C e ■ if = f e 
1 d t v° + C° ■ v° = r ° 



with 



q e - MJ° • D x u° 



M™ ■ D y v° 
M° e • D y iT 



(15) 



where C e , C°, f e , and r° are time- independent. Moreover, since the matrix C = 
diag(ci, . . . ,cj\f) is diagonal, the equations in (15) decouple into scalar equations of the form 

d T u k (x,y,T) + c k (x,y)u k (x,y,T) = f k (x,y) (16) 

that need to be solved from r = t until r = t + At. In ( Jl6| ), we have c k {x, y) = c k (x, y,t + ^At) 
and f k (x,y) = r k (x,y,t + ^At). The exact solution of (16) is 

u k (x, y, t + At) = exp(-c fc (x, y)At)u k (x, y, t) - g^r^y (1 - exp(-c fc (x, y)At))f k (x, y) (17) 

= u k (x, y, t) + At (r k (x, y) - c k (x, y)u k (x, y, t)) E(-c k (x, y)At) , (18) 

where E(c) = ""'f 1 . Note that, given a robust implementation of this function E (see 
Sect. [5]), the representation (18) has an important advantage over formula (17), namely it 
is defined in voids, where c k (x,y) = 0. Using formula (18), we therefore define the even 
evolution operator 

ei = [si + ^(-c^) r ■ e ' i 



(19) 



as the one that advances the even components according to (18), and the odd evolution 
operator 





CO 

°t+At,t 



+ AtE(-C°At) 



(20) 



as the one that advances the odd components according to (18). A full solution step from t 
to t + At is then achieved via four sub-half-steps 



St- 



•At,t 



00 

o. . 1 



Sf 



s?. 



(21) 



Note that in general Sf i A o S° 



t+^Att ^ ^t+Atf> as one can see f rom t ne solution formula 



(18). 



4.4. Accuracy. Due to the proper setup of the solution components on staggered grids, all 
spatial differential operators are approximated in a central fashion, and thus with second order 
accuracy, i.e., the spatial truncation error is 0(h 2 ), where h = max{Ax, Ay}. Moreover, due 
to the symmetry in the arrangement of the sub-half-steps in (21), the splitting error in a 
single time step is 0(At 3 ), and henceforth the global temporal truncation error due to the 
fractional steps is 0(At 2 ) [22\. This second order accuracy in time is preserved because all 
temporal evaluations are done in a symmetric fashion (at the half-step time t + \ At). The 
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overall second order accuracy (i.e., the truncation error is 0(At 2 ) + 0(h 2 )) is confirmed by 
the numerical results in Sect. 16.11 



4.5. Stability. Here, we show the L 2 stability at the heart of the approach, namely the 
combination of staggered spatial grids with the temporal splitting ( pTj ). We conduct the 
analysis in a simplified case, namely: the Pi system (with advection coefficients set to 1), in 
one space dimension, with constant coefficients, and without a source (i.e., q = 0). In this case 



our solution vector consists of two scalar fields u = \i 

{dtu e + d x u° + c e u e 
d t u° + d x u c + c°u° 



u \ 

- 
= . 



which satisfy the equations 



We now conduct a von Neumann stability analysis of the numerical scheme presented above. 
To that end, let the numerical solution at time t be represented in terms of its Fourier 
coefficients 

u e (x,t) = ^2a e k (t)e ikx and u°(x, t) = ^ a° k {t)e ikx . 

k k 

Note that this particular form of the Fourier representation holds for a periodic domain of 
length 2ir. However, the arguments below transfer to any periodic domain, as well as to an 
infinite domain (Cauchy problem, for which the sums would turn into integrals). 



A step of the numerical scheme (21 ) in general affects all the Fourier coefficients a k and a k . 
However, due to the linearity of the update rule, no mixing occurs between Fourier modes 
of different frequency, i.e., the coefficients a e k (t + At) and a k (t + At) depend only on the 
coefficients a k (t) and a k (t) for this particular k. Consequently, it suffices to investigate the 
growth factor for basic wave solutions 



ikx 



u e (x, t) = a e k (t)e tkx and u°(x, t) = a° k {t)e 

Let u e be defined on the grid hj, j € Z and u° be defined on the grid h (j + |), j G 
the (staggered) grid solution values are 



(22) 
7 j. Then 



u c {hj,t) = a%(t)e ikhj and u°(h{j + 
and consequently the staggered grid derivatives are 



I). 



t) = a° k {ty kh{ i+k) 



and 



(D x u e )(h(j + ±),t) = a e k (t) 



(D x u°)(hj,t) = a° k (t) 



ikh (j+l) 



ikhj 



h 



2^sm(kh/2)e lkh(j+ ^a c k (t) 



j,kh(j+\) _ e ikh(j-i) 
h 



2^sm(kh/2)e ikhj a° k (t) 



Hence, for the basic waves (22) on the staggered grids, the even half-step solution operator 



( 19 ) updates the Fourier coefficients by the linear operation 

(t + I At) 



4 





r 









i 
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exp(— hfAt) and f e 



where d e 

solution operator (20) acts on the Fourier coefficients as 

11 (< + \ Ai) 



i^-E(-\c e At)sm(kh/2). Similarly, the odd half-step 
(*), 



" 1 o" 




'4 


f° d °_ 




<_ 



where d° = exp(— |c°Ai) and f° = —i^E(—^c°At)sm(kh/2). Consequently, a full solution 



time step ( [21] ) acts on the Fourier coefficients as 
(t + At) = G°i At ■ Gl Af ■ G\ At ■ G°i At 



4 



ai 



(d e y + (d e + i)/ e /° 



d°(d e + l)/ e 



(d° + (d e y)f° + (d e + i)f e (f°y (dry + d°(d e + i)/ e /° 





~ a t~ 







(t) 



=G 

The eigenvalues of the growth factor matrix G are 

Al,2 



9 



CrlO\2 



(d c d' 



where g = \ {{d e ) 2 + (d°) 2 + (d e + l)(d° + l)/ e /°) is half the trace of G. Since d e and d° are 
real, and f e and /° are purely imaginary, the half-trace g is real. 



Theorem 3. If At < h, then | Ai 2 1 < 1; *- e v ^ e iiwie stepping (21) yields a stable scheme. 



Proof. First, one can observe that if |g| < d°d°, then both eigenvalues Ai^ are complex, with 
real part g, and imaginary part ±\/ (d c d°) 2 — g 2 . Therefore, their magnitude satisfies 

|Ai, 2 | 2 = g 2 + {(d c d°) 2 - g 2 ) = (d e d°) 2 < 1 , 
i.e., the stability condition is satisfied. 

Second, if g > d e d° > 0, then by directly estimating g we have 

g 2 - (d e d°) 2 = (g- d c d°)(g + d e d°) < \{{d c ) 2 + (d°) 2 - 2d c d°) |((d c ) 2 + (d°) 2 + 2d c d°) 
= l(d e -d ) 2 (d c + d°) 2 , 

and therefore 

A 2 = 9 + yV - (d c d°) 2 < \{{d e ) 2 + (d°) 2 + \d e - d°\ (d e + d°)) = max{(d e ) 2 , (d°) 2 } < 1 . 

Moreover, since g > 0, we have that Ai > g — g = 0, hence the stability condition is satisfied. 

Finally, we need to consider the case g < —d e d° < 0. Clearly, A 2 < g + \g\ = 0, thus 
stability is satisfied if Ai > —1, which is equivalent to the condition 1 + (d e d°) 2 + 2g > 0. The 
half-trace g depends on the wave number k. It achieves its smallest value if sin(fc/i/2) 2 = 1, 
thus we can estimate (using the CFL condition At/h < 1) 

2g > {d c ) 2 + {d°) 2 - (1 + d c )E{-\c c At) (1 + d°)E{ 



c°At) 



■o\2 



(d e f + (d°) 



io\2 



1 - (d e ) 2 1 - (d°) 



\<?At 



\c°At 



exp(-c e At) + exp(-c°At) - 4 
2 



1 - exp(-c e At) 1 - exp(-c°At) 



c e At c°At 
c e AtE{-c e At) + c°AtE{-c°At) + 4E(-c e At)E(-c°At)) 
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and therefore (using the notation t e = c e At and t° = c°At) we obtain 

1 + (d c d°) 2 + 2g>4- 2fE(-f) - 2t°E(-t°) + (ft° - 4)E(-f)E(-t°) > , 

where this last inequality can be verified by expanding the power series of the function E. This 
shows that the presented scheme is conditionally stable, with the CFL condition At < h. □ 

In analogy with similar types of problems and numerical approaches, one can expect that 
these I? stability results carry over to the general case, albeit with the following modifications: 

• By Duhamel's principle, problems with source terms can be interpreted as a super- 
position of many homogeneous problems with new initial values. Hence, the stability 
results carry over to equations with sources. 

• For systems with more than two components that pre-multiply the spatial derivative 
by a matrix M, the maximum admissible time step has to be scaled by the inverse of 
the largest (in magnitude) eigenvalue of M. 

• In two space dimensions, an extra factor of h in front of the time step will account 
for the presence of growth rates in each of the two spatial dimensions. 

• The case of variable coefficients is not covered by von-Neumann analysis. However, 
as the proof of stability shows, larger decay terms relax the requirements for stability. 
Hence, it is plausible that the variable coefficient case does not exhibit worse stability 
properties than a constant coefficient case with the same minimal decay rates. 

The numerical results in Sect. [6] confirm the stability of the method in the general case. 



4.6. Important Advantages of the Methodology. Besides its aforementioned second 
order accuracy (in space and time), an important advantage of the presented approach is 
its structural simplicity and regularity: in each time step and at every grid point, exactly 
the same operations are performed (albeit with different coefficients), thus giving rise to 
an efficient vectorization (see Sect. [5]), and possibly parallelization (see Outlook). Another 
advantage of the approach is its stable treatment of large decay coefficients (i.e., problems 
with large absorption and/or scattering): due to the exact solution (|18|) of the sub-step ODE 



( 16 ), large values of the decay coefficients do not impose restrictions on the time step. This 
is in contrast to explicit time-stepping rules (such as Runge-Kutta methods), which would 
incur severe time step restriction for such stiff problems. 



4.7. Limitations of the Methodology. The simplicity of the approach incurs a few fun- 
damental limitations, and the user of StaRMAP should be aware of these limitations when 
using the code. Most prominently, the use of the P/v closure (or SPn) gives rise to the Gibbs 
phenomenon (see Sect. [TJ, i.e., the exact solution of the moment system may develop oscil- 
lations that are absent in the solution of the original radiative transfer equation ([T|. As one 
consequence, the vector of moments may not be realizable, i.e., the moments cannot stem 
from a non-negative density. This phenomenon is well-known to users of spherical harmonics 



moment equations, and it is particularly demonstrated in the line source test case in Sect. 6.3 



There is also a Gibbs phenomenon due to having a linear second order numerical scheme, 
which due to Godunov's limitation theorem [6] cannot be monotone, and thus spurious oscil- 
lations tend to occur near strong gradients of the solution. These overshoots could be avoided 
through limiters. However, the addition of limiters would go at the expense of the structural 
simplicity and efficiency of the method, and they are therefore not included in StaRMAP. 
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Another limitation of the approach is that its simplicity stems from the regularity in the 
geometry. Hence, a generalization of the method to irregular or locally refined meshes is 
not possible without sacrificing some fundamental structural advantages. Similarly, since the 
method's overall second order accuracy is based on exploiting local symmetries, a generaliza- 
tion to higher orders is impossible without introducing major modifications. 

Finally, the temporal splitting could generate spurious oscillations in the case of strong 
spatial gradients in the material parameters. While the exact treatment of the sub-step ODE 
( |16| ) successfully deals with uniformly large decay rates, extreme gradients in the absorption 
and/or scattering coefficients could trigger O(l) spurious waves, if the time step is not suitably 
reduced. Again, it is important that the user of StaRMAP be aware of this possibility when 
attempting to solve problems with large material parameters that in addition exhibit strong 
gradients or jumps. 



5. Implementation in Matlab 

At the core of the Matlab implementation of the method developed in Sect. [4] is the 
fact that the differential and evolution operators perform the same type of operation at each 
grid point. We therefore store each field quantity and each component of the solution vector 
as a two-dimensional array (which is identical to a matrix in Matlab), in which the first 
component represents the x-index and the second component the y-index of the respective 
grid. Hence, in line with the staggered grids ^ and (11), every component in the set In is a 



matrix of size n x x n y , and the components in I^i, I\2, and I22 are of sizes (n x J rl—p x j ~ '<<?/, 
n x x (n y + 1 — p y ), and (n x + 1 — p x ) x (n y + 1 — p y ), respectively, where p x and p y are the 
periodicity flags defined in ( Jl~2| ). 

5.1. Spatial Derivatives. In the StaRMAP code, the moments are assembled in a cell struc- 
ture U, where the j component is U{j}. Note that components that are defined on different 



grids may be of different sizes. The half-grid central difference formulas (13) applied to the 
component U{j} are therefore computed as 

dxTHj} = diff(U{j}(extendx{l},:),[],l)/li(l); 
dylHj} = diff(U{j}(:,extendy{l»,[],2)/h(2); 

where the vector h contains the grid spacing Ax and Ay as components. As described in 
Sect. 4.1, if the component U{j} is defined on the grid Gki, then dxU{j} is defined on Gz-k,e, 



and dyU{ j} is defined on Gk.3~£- Since the dif f function reduces one of the array's dimension 
by 1, the array U{j} may need to be extended, before finite differencing is applied. This is 
encoded in the index vectors extendxjl}, extendx{2}, extendyjl}, and extendy{2}. The 
definition of these four index vectors is also the only place in the code where the boundary 
conditions enter. 

5.2. Update of Even and Odd Components. Each time step consists of four half steps. 



The even solution operator (19) updates the components on the even grids G\\ and G22, by 



using spatial derivatives of quantities on the odd grids G21 and G21] and vice versa for the 



odd solution operator (20). One step with the even (odd) solution operator consists therefore 
of two parts: first, spatial differences of the odd (even) components are computed (see above); 
second, the even (odd) components are updated, via the code segment 
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W = -sumcell([dxU(Ix{j}) ,dyU(Iy{j})] , . . . 

[par.Mx(j ,Ix{j}) ,par .My(j , Iy{j})] ) ; 
U{j} = U{j}+dt/2*(W+Q{j>-. . . 

sT{gtx(j) ,gty(j) ,par.mom_order(j)}.*U{j}) . *. . . 

ET{gtx ( j ) , gty ( j ) , par . mom_order ( j ) > ; 

Here the array W contains all D x and D y finite differences that influence a certain component, 
weighted by the negative moment matrices M x and M y . The call of the function sumcell 
achieves this task efficiently by only considering derivatives that correspond to actual nonzero 
entries in the matrices M x and M y . Since the moment matrices derived in Sect. [3] have an 
O(l) number of nonzero entries per row, this guarantees that the computational effort of the 
code grows linearly with the number of moments. 

The update rule for the j th moment is an immediate implementation of the update for- 
mula (18). The array Q{j} represents the j th moment of the source term, and the arrays 



sTjgtx ( j ) , gty ( j ) , par . mom_order ( j ) } and ET{gtx ( j ) , gty ( j ) , par . mom_order ( j ) } encode 
the decay quantity Cj(x,y) and E(—Cj(x,y)^), respectively. Here, sT stands for Sj, and ET 
denotes the function E being applied to Sj. There is a small modification to this update rule 
for j = 1: since the zeroth moment is unaffected by scattering, the update is conducted with 
the fields sA (which stands for E a ) and EA (which is the function E applied to 

Structurally, the decay terms sT and ET are arranged in three dimensional cell arrays. 
The first two components encode the grid index in each coordinate direction (1 for odd, and 
2 for even). The third component encodes the moment order index £, given in the index 
vector par .mom.order, which encodes a mapping from the system component index to the 
moment order t. Since different moments of the same moment order t possess the same decay 
parameters, this structure guarantees that the number of decay quantities grows only linearly 
in the moment order N. 



5.3. expmldiv. The function E(c) = exp ^~ 1 , used in the solution formula (18), requires a 
careful implementation. For values of c away from zero, this formula can be implemented as 
given. However, for |c| <C 1 the difference and division may lead to severely amplified round- 



off errors, which to leading order equal -, where S is the machine accuracy. In the StaRMAP 



c 



code (function expmldiv), the exponential formula is therefore replaced by the Taylor series 
approximation E(c) = 1 — \c-\- |c 2 for |c| < 2 x 1CP 4 . While this series approximation does 
not yield any amplified round-off errors, it possesses an approximation error that to leading 
order equals 53 c 3 . For double-precision arithmetics, the total errors incurred with these two 
evaluation formulas match for |c| ~ 2 x 10 -4 , at a value of less than 10 -12 . 



5.4. Evaluation of Material Parameters and Source. In order for the aforementioned 
update step to have the decay quantities sA and sT, as well as the source Q available, these 
quantities must be evaluated before. While the StaRMAP code can treat rather general 
problem setups (anisotropic and/or time-dependent parameters), an important focus lies on 
its efficiency in special situations, and a significant part of the code is devoted to properly 
addressing this demand, as follows: 

• Scattering, absorption, and the source term are each evaluated once in the beginning 
of each time step, at time t + ^Ai. However, if any of these quantities is actually 
time-independent, then its evaluation occurs in the first time step only. 
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• The scattering is divided into its isotropic component (encoded in par . sigma.sO) 
and the anisotropic component (encoded in par . sigma.sm), where m stands for the 
moment order (see above). For a problem with isotropic scattering, one simply does 
not provide the function par . sigma_sm. In turn, if this function is provided, then it 
is evaluated for all moment orders 1, . . . ,N. Note that the function par. sigma.sm is 
never evaluated for m = 0; scattering never influences the zeroth moment. 

• The source function par. source allows for a component for each moment. However, 
since in many applications sources are acting on only a few moments, one can encode 
the indices of these moments in the vector par . source_ind. Specifically, an isotropic 
source would apply only to the zeroth moment, in which case par . source_ind = 1. 

5.5. Time Stepping. Each time step consists of three parts: (i) the evaluation of mate- 
rial parameters and sources, if applicable; (ii) the application of the four half-step solution 



operators, given in (21); and (iii) the plotting of the solution, if applicable. 

In the application of the half-step operators, the even solution operator S* i ^ f is applied 

twice in succession. While one cannot replace these two half-steps by one full step, the finite 
differences need to be evaluated only once. This is done in the StaRMAP code, thus reducing 
the number of operations. Note that in the special case of a fully time-independent problem, 
one could apply the same trick to the odd solution operator (combined over two successive 
time steps). Since this case is quite rare, this trick is not implemented. 

The (zeroth moment of the) solution is plotted at all times given in the vector par . t_plot. 
At the end of each time step, for all times in par .t_plot that lie in the interval (t,t + At], the 
solution is computed via linear interpolation in time u(t p \ t) = (1 — X)u(t) + Xu(t + At) , where 
A = - £ ^ r 1 . This preserves the second order accuracy of the solution. With this approach, 
the computation is completely unaffected by the output of the solution. 

5.6. Initialization. A significant part of the StaRMAP code is devoted to the initialization of 
the data structures. While many of these steps are technical (mostly to ensure that Matlab 
handles the memory in an efficient fashion), some steps are of conceptual importance: 

• The staggered grid index sets ell, c22, c21, and cl2 (representing the respective 
sets In, I21, I12, and I22 are created automatically from the structure of the moment 
matrices par . Mx and par . My. This is done by first placing the first solution component 
in ell, and then keeping on placing components in appropriate sets, whenever a 
nonzero entry in the matrices par . Mx and par . My requires this, until all components 
are distributed into the index sets. 

• The maximal admissible time step is determined by computing the largest (in mag- 
nitude) eigenvalue of the moment matrix M x via the expression 

abs(eigs (par .Mx, 1 , ' lm' ) ) 

Note that it is assumed that the moment system preserves rotational symmetry, and 
thus the eigenvalues of M x are identical to the eigenvalues of n x M x + n y M y for all 
n 2 x + n 2 y = 1. 

• The code segment 

extendx = { [1 :par .be (1) , 1 :nl (1) ,par .be (1) * (nl (1) -1)+1] , . . . 

[n2(l)*(l:l-par.bc(l)) ,l:n2(l)]}; 
extendy = { [1 :par .be (2) , 1 :nl (2) ,par .be (2) * (nl (2) -1)+1] , . . . 

[n2(2)*(l:l-par.bc(2)) , 1 :n2(2)] }; 
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creates the index vectors that encode the boundary conditions, as described above. 
• In the case of isotropic scattering, the vector par .mom.order is modified to point 
at the first moment order only, resulting in the three-dimensional cell array of the 
material parameters to be of length 1 in the third component. 

5.7. Definition of a Test Case. The StaRMAP project is designed so that in all "usual" 
cases, the file starmap.solver .m does not require any modification. Each test case is encoded 
in an example file (such as starmap_ex_checkerboard.m), which defines the problem and 
then calls the solver file. The communication is achieved via the Matlab struct par, which 
is created in the example file and then passed to the solver file. This problem parameter 
can contain the problem name (par. name), the moment matrices (par.Mx and par. My), the 
moment order index vector (par .mom_order), the coordinates of the computational domain 
(par . ax) , the numbers of grid cells in each coordinate direction (par . n) , the type of boundary 
conditions (par. be), the times of plotting (par . t_plot), the final time (par .tf inal), the 
CFL number (par.cn), as well as function handles for absorption (par . sigma_a), scattering 
(par . sigma_sO and par . sigma_sm), the source (par . source), initial conditions (par.ic), 
and a problem-specific plotting routine (par . output). Any of these parameters that is not 
provided is set to a default value by the solver file. 

A typical StaRMAP example file consists of three parts: (i) the definition of the problem 
parameters that differ from the default, including function handles; (ii) the creation of the 
moment matrices (by calling starmap_closure_pn.m for P/v and starmap_closure_spn.m for 
SPn) and the call of the solver; and (iii) the definition of the problem-specific functions. 



6. Numerical Results 

6.1. Verification via Manufactured Solutions. We use the method of manufactured so- 
lutions [19] to verify our implementation and to validate the accuracy predictions about the 



numerical approach made in Sect. 4.4 To that end, a routine is implemented that uses Mat- 
lab's symbolic toolbox to automatically compute the source vector q(x, y, t) that generates a 
prescribed solution under prescribed material coefficients. We select a smooth, but spatially 
and temporally varying, solution and smooth, but non-constant, coefficients. Moreover, the 
solution and the coefficients (and consequently the sources) are periodic on the rectangular 
computational domain. 

Figure [2] displays the spatial L 1 , L 2 and L°° errors in the scalar flux R® at the final 
time as a function of the spatial resolution of the discretization. It shows the expected 
second order convergence. In this case, we consider the P3 model for t £ [0, 0.5] on the 
domain (x,y) G [0,1] x [0,1]. We use a time- and space-dependent absorption coefficient 
T, a (t,x) = t cos(27TX2), and anisotropic Henyey-Greenstein scattering [5] with T, s £ = 0.9 e . The 
manufactured solution is 

Ro(t,x) = e _t sin 2 (2vrxi) 
with all other moments being zero. 

In the StaRMAP project, the "manufacturing" of a solution, i.e., the computation of a 
source term that generates a prescribed solution, in done via the file starmap_createjnms .m. 
This m-file generates the actual example file starmap_exjnms_auto .m, which then can be 
executed to produce the plot shown in Fig. [2} 
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Manufactured Solution Test with P3 




2 4 6 



Figure 3. Checkerboard problem, material coefficients: isotropic source (white), purely scattering 
E s n = 1 — T t (orange and white), purely absorbing T a = 10 = T t . 

6.2. Checkerboard. We consider the checkerboard problem described in [I]. The compu- 
tational domain is a square of size [0, 7] x [0, 7] where the majority of the region is purely 
scattering. In the middle of the lattice system, in the square [3, 4] x [3, 4], an isotropic source 
q = 1 is continuously generating particles. Additionally, there are eleven small squares of side 
length 1 of purely absorbing spots in which T a = 10 = Sj. In the rest of the domain, T a = 0, 
Tt — 1- Figure [3] illustrates the problem setting more precisely. 
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Checkerboard Test with P3 at t = 3.20 



Checkerboard Test with P5 at t = 3.20 




1 2 3 4 5 6 
Checkerboard Test with P15 at t = 3.20 



1 2 3 4 5 6 
Checkerboard Test with P39 at t = 3.20 
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Figure 4. Checkerboard problem: Scalar flux at t — 3.2 for several moment approximations (P3, 
P5, P15, and -P39), computed with 250 x 250 grid points. The values are plotted in a logarithmic scale 
and limited to seven orders of magnitude. 



In the StaRMAP project, the checkerboard test case is implemented in the example file 
starmap_ex_checkerboard.m. Extrapolation boundary condition^] are enforced. At the ini- 
tial time t = 0, all quantities are zero. The spatial grid is 250 x 250 for all cases. 

We compare the scalar flux using different orders of P/v approximations, with the scalar 
flux Rq at t = 3.2 shown in Fig. g The case N = 3 corresponds to 10 angular degrees of 
freedom, while the case A^ = 39 possesses 820 angular degrees of freedom. For increasing N, 
several well-known properties of the P/v approximation can be observed. First, the maximum 
propagation speed of information approaches the correct limit of 1 which can be seen at the 



Most frequently, this test is equipped with vacuum boundary conditions, which are not implemented in 
the StaRMAP routines. The easiest way to "emulate" vacuum in StaRMAP for this test would be to extend 
the computational domain. 
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Line Source Test with P19 at t = 0.50 
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Line Source Test with P39 at t = 0.50 
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Figure 5. Results of the line source problem at t = 0.5, computed with Pig (top) and P^g (bottom). 
The left sub-figures show the 2d solution profile. The right sub-figures show the solution on a cut 
along the positive radial axis (red), together with Ganapol's semi-analytic solution (blue). 



upper front. Second, the Gibbs phenomena in the solution (left, right, bottom) are diminished. 
The P39 solution can be seen as a fully converged transport solution. 

To get an impression about the computational cost and efficiency of the StaRMAP solver: 
the computation of a P5 solution on a 100 x 100 spatial grid corresponds to 18.3 million 
degrees of freedom in space, angle, and time, and it takes about 1.5 seconds to compute on a 
2008 Lenovo W500 laptop. 



6.3. Line Source. The line source problem has been investigated for the Pn equations in 
[lj. Essentially one is trying to compute the Green's function for an initial isotropic Dirac 
mass at the origin, i.e. 

ip(t = 0,x,U) = —Six) . 
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One can also imagine an infinite line emitting particles isotropically (hence the name). There 
exists a semi-analytical solution to the full transport equation ([I]) for this problem due to 
Ganapol et al. [5] . The exact solution consists of a circular front moving away from the origin 
as well as a tail of particles which have been scattered or not emitted perpendicularly from 
the line. The front does not consist of Diracs, but has a singularity at the edge. 

We choose £ a = and S s = so that £ s o = 1. In our numerical experiments, the initial 
condition is realized as a narrow Gaussian in space 

R = JZZ ex P 



with a = 3.2 x 10 -4 . This way, the integral of Ft® over the whole space is one, so it can be 
directly compared to the scalar flux of the benchmark solution. The computational domain 
is [-0.6,0.6] x [-0.6,0.6]. 

In the StaRMAP project, the line source test case is implemented in the example file 
starmap_ex_linesource .m. This m-file also approximates Ganapol's semi-analytic solution 
with an accuracy that makes the resulting function indistinguishable from the exact solution 
when plotted. The computational results for Pig and P39, both computed on a 150 x 150 
spatial grid, are shown in Figure [5] These results are in line with the observations from [1] 
and the previous test case. While in theory the P/v equations have a rotationally symmetric 
solution, it has been observed in pQ that a higher-order scheme is necessary to observe the 
rotational invariance numerically. In our case, both the P19 and the P39 solutions exhibit this 
behavior. For the lower-order P19 model, the Gibbs phenomenon is very clearly visible. The 
P39 solution, on the other hand, shows only very tiny oscillations. Of course, it cannot capture 
the singular behavior of the solution, but it does a good job at capturing the radiation front. 
These two examples are shown because the StaRMAP code solves them within a few minutes. 
However, if one is willing to accept longer computation times, then suitable grid refinement 
and a more narrowly focused initial condition would yield an even better approximation. 



6.4. Beam in Void and Medium. In this problem we study the advection of a packet of 
particles that are emanating from a source (essentially compactly supported in space) and all 
move in the same direction in empty space, before they hit a material. This test investigates 
how well the method performs in a void, at material interfaces and with anisotropic scattering. 

The source can be written as 

q(t, x, f2) = q(x)5(£l — f^o) , 

where q is a spatial distribution and £Iq is a direction in the plane. Here we choose f^o to 
have an angle of tt/6 with respect to the x-axis and q is the same spatial Gaussian as for 
the line source problem. The proper solution of this beam problem involves the challenge of 
calculating the correct moments of the source for all components of the Pat solution vector. 
We do so by following the definitions and transformations given in Sect. [3j 

The domain for this test is [—0.6, 0.6] x [—0.6, 0.6] with extrapolation boundary conditions 
on all sides. The region with x < 0.3 contains no material, i.e. £j = 0, whereas the region 
x > 0.3 is a material with no absorption (£ a = 0) but anisotropic Henyey-Greenstein [8] 
scattering T, s £ = W0g e with g = 0.85. 

In the StaRMAP project, this test case is realized similarly to the manufactured solution 
example, namely the file starmap_create_beam.m performs the computation of the beam 
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Beam Test with P9 at t = 0.30 Beam Test with P9 at t = 0.60 




Beam Test with P39 at t = 0.30 Beam Test with P39 at t = 0.60 




Figure 6. Beam problem, computed with Pg (top row) and P39 (bottom row). Shown is the scalar 
flux i?° at times t = 0.3 (left column) and t — 0.6 (right column). A spatially Gaussian source 
propagates into the direction with angle ir/6 with respect to x-axis, and hits a medium at x = 0.3. 

initial conditions, and then generates an example file starmap_ex_beam_auto .m that can then 
be executed to produce the plots shown in Fig. [6j 

Figure [6] shows the scalar flux R® at two times (t £ {0.3,0.6}) for a beam computation 
on a spatial grid of 150 x 150 cells, using two moment models: Pg (top two figures) and 
P39 (bottom two figures). The full transport equation would initially (i.e., before the 
beam hits the medium) just advect the source with unit speed into direction Oo- Due to the 
truncation in the P/v approximation, the beam smears out, and Gibbs artifacts appear during 
the evolution. These are clearly visible in the Pg solution (top left figure), in which the beam 
spreads out and oscillations occur. In contrast, the P3g solution at t = 0.3 (bottom left figure) 
is very close to a fully converged transport solution; no truncation artifacts are visible. 

At the material interface, particles start to scatter strongly, albeit in a forward direction. 
As a consequence, the propagation of the beam effectively slows down and it is smeared out. 
The P3g solution at t = 0.6 (bottom right figure) shows the correct beam profile (including a 
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reflected beam at the interface between medium and void), albeit a small amount of Gibbs 
artifacts that are visible. The corresponding Pg solution (top right figure) is of significantly 
lower accuracy: while the general direction of the radiative flux is represented correctly, the 
precise shape of the radiation profile (e.g., beam thickness, reflections, maxima) is not well- 
captured. 

When considering this test, it must be remarked that beams are not the type of phenomena 
that one would usually consider using moment methods for. However, as the results in Fig. [6] 
show, and as one can easily verify with the StaRMAP project files, with the moment order 
N chosen suitably large, one can in fact compute beams with a high accuracy. Therefore 
challenging problems, such as beams moving from a void into a highly scattering medium, 
can be computed accurately with StaRMAP. 

6.5. Further Test Cases. The StaRMAP code has also been applied to several other exam- 
ples. Among them are: 

• The "boxes" test, implemented in the example file starmap_ex_boxes .m, is a test case 
that was proposed by McClarren [12] to demonstrate the equivalence of the SPjsr and 
the P/v equations under certain assumptions. 

• The "control rod" test was devised by Olbrant |17] as an example for time-dependent 
material coefficients. The setup can be interpreted as a simple model for the evolution 
of the neutron density in a nuclear reactor, when an absorbing rod is moved into and 
out of the domain. In the StaRMAP project, it is implemented in the example file 
starmap_ex_controlrod . m. 

Further details and numerical results for both tests are given in [17]. 



7. Conclusions and Outlook 

The presented numerical approach to solve the spherical harmonics P/v equations (and 
variants thereof) of radiative transfer has been demonstrated to be highly flexible, and to 
be applicable to a wide variety of radiative transfer problems in two space dimensions. The 
method is implemented in the Matlab project StaRMAP, that can be downloaded |21j . and 
all examples presented in this paper can be reproduced by the user. 

As shown in this paper, the staggered grid finite difference approach employed in StaRMAP 
is second order accurate, and stable independent of the magnitude of the material scattering 
and absorption coefficients (with the CFL condition that is commonly incurred for advection 
problems). Further advantages and drawbacks of the method have been discussed in detail. 

The availability of the complete Matlab code serves another purpose besides the repro- 
ducibility of the results presented in this paper: we encourage users to create their own 
StaRMAP example files for test cases that arise in their own work. The modular structure of 
the code facilitates this goal: in general, users should not be required to modify the solver 
file starmap_solver .m, but instead they can simply build on the existing example files. 

While the current StaRMAP code is quite general in that it allows for arbitrary moment 
order and for various types of moment systems, it has a number of crucial limitations. Some 
cannot be overcome easily, such as the occurrence of Gibbs phenomena or the regularity of 
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the grid structures (see Sect. 4.7 for more details). In turn, some limitations can (and shall) 



be addressed in the future, as follows. 

In principle, the presented methodology applies to the full 3D case as well. We have 
decided against generalizing the current StaRMAP code to 3D, mainly for the simple reason of 
computational cost and memory requirements. While 2D problems can be solved in Matlab 
with a good accuracy, in reasonable compute times, and with a reasonable amount of memory 
consumption, we do not perceive things as quite mature enough yet for a 3D Matlab version 
of StaRMAP. 

The StaRMAP code can be extended to other moment models that exhibit a similar coupling 
structure between the moments as the P/v and the SPn equations. Specifically, the filtered P/v 
(-F-P/v) method [H] falls into this category, as do the diffusion-corrected Pn (Dn) equations, 
which contain an additional diffusion term in the highest order moments. These last types of 
models can be rationalized via asymptotic analysis [20] or by the Mori-Zwanzig formalism of 
statistical mechanics [3]. The splitting procedure employed in the StaRMAP code allows for 
the addition of a diffusion step applied to the moments of the highest order. Finally, it must 
be verified that the method presented here satisfies an important property for numerical 
approaches for radiative transfer, namely whether it is asymptotic-preserving (AP) [9]. It 
can be quite easily shown that the scheme proposed in this paper is at least formally AP. 
Moreover, as has been suggested in [10], because of the staggered grid we obtain a compact 
stencil in the diffusive limit. In addition, due to the explicit integration in time, the time step 
does not depend on the possibly stiff right hand side. The AP property will be investigated 
in detail in a companion paper [3]. 
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